﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary
{
    /// <summary>
    /// This class provides the approximation of the Jacobian matrix at a specified position.
    /// </summary>
    [Serializable]
    public class JacobianMatrixApproximator
    {
        /// <summary>
        /// This field holds the system of equation.
        /// </summary>
        private NonLinearEquationSystem equationSystem;

        /// <summary>
        /// Initializes a new instance of the <see cref="JacobianMatrixApproximator"/> class.
        /// </summary>
        /// <param name="equationSystem">The nonlinear system of equation.</param>
        public JacobianMatrixApproximator(NonLinearEquationSystem equationSystem)
        {
            if (equationSystem == (NonLinearEquationSystem)null)
                throw new ArgumentNullException("equationSystem");

            this.equationSystem = equationSystem;
        }

        /// <summary>
        /// Gets or sets the nonlinear system of equation.
        /// </summary>
        /// <value>The nonlinear system of equation.</value>
        public NonLinearEquationSystem EquationSystem
        {
            get { return this.equationSystem; }
            set { this.equationSystem = value; }
        }

        /// <summary>
        /// Approximates the Jacobian matrix at a specified position.
        /// </summary>
        /// <param name="positionVector">The specified position for which the Jacobian matrix should be approximated.</param>
        /// <param name="precision">The precision of the approximation.</param>
        /// <returns>The approximated Jacobian matrix at the specified position.</returns>
        public Matrix ApproximateJacobianMatrix(GeneralVector positionVector, double precision)
        {
            double tmp,delta;
            int dimension = this.equationSystem.EquationSystem.Length;
            GeneralVector x = positionVector.Copy();
            Matrix result = new Matrix(dimension, x.Count);
            GeneralVector fv = this.equationSystem.SolveAt(x);

            for (int i = 0; i < x.Count; i++)
            {
                tmp = x[i];
               
                if (tmp > 1.0)
                {
                    delta = precision*tmp;
                }
                else
                {
                    delta = precision;
                }

                x[i] = tmp + delta;       
                delta = x[i] - tmp;      

                GeneralVector p = this.equationSystem.SolveAt(x);
              
                x[i] = tmp;

                for (int j = 0; j < dimension; j++)
                {
                    result.SetValueAtPosition(j,i,  (p[j] - fv[j]) / delta);
                }
            }

            return result;
        }
    }
}
